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. Abstract 

In a variety of applications, most notably microfluidic design, slip-based boundary conditions have been sought 
CN| ' to characterize fluid flow over patterned surfaces. We focus on laminar shear flows over surfaces with periodic 
height fluctuations and/or fluctuating Navier scalar slip properties. We derive a general formula for the "effective 
^ . slip" , which describes equivalent fluid motion at the mean surface as depicted by the linear velocity profile that 
1^ . arises far from it. We show that the slip and the applied stress are related linearly through a tensorial mobility 
^ I matrix, and the method of domain perturbation is then used to derive an approximate formula for the mobility 
^ I I I law directly in terms of surface properties. The specific accuracy of the approximation is detailed, and the 
C/3 ' mobility relation is then utilized to address several questions, such as the determination of optimal surface 
• i-H ' shapes and the effect of random surface fluctuations on fluid slip. 
>^ 
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I 1. Introduction 
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With recent progress in the study and fabrication of mir ofluidic devic e s, new interest has arisen in determinin g 



J> ' generalized forms of hydrodynamic boundary conditions (jStone et al.l . l2004l : iBazant and Vinogradoval . l2008l ) 

OO ■ 



Advances in lithography to pattern substrates at the micrometer and nanometer length scales have raised 

CO ' several questions in the modeling of fluid motions over these surfaces. For example, rather than trying to solve 

equations of motion for the flow at the scale of the individual corrugations of the pattern, it is appropriate to 

consider the bulk fluid motion (on length scales much larger than the pattern wavelength) by utilizing effective 

On '. boundary conditions that characterize the flow at the surface. For fluid being sheared horizontally over a 

textured surface, a basic aim is to derive a local boundary condition that can be applied along the smooth 

mean surface, which mimics the effects of the actual condition along the true surface. These effective conditions 

' can be used in place of the no-slip condition to solve for macro-scale flow without the tedium of enforcing a 
?— ( ' 

' boundary condition on a rough boundary geometry. 

A standard phenomenological approach for flow over a patterned surface is to assume a Navier slip boundary 
condition 

u^=U-u = 6^ (1) 
on 

which relates the fluid velocity u at the surface, the velocity of the surface U, and the shear strain rate normal 
to the mean surface, du/dn, via the slip-length b. This approach has been studied extensively in experiments, 
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theore t ical calculations, and simul at ions (for recent discussio n and results see I Vinogradoval (jl999l ): iLauga et al 
(|2007h : lBocauet and BarratI pOOTT ) : loavis and Laueal |2009[ )). 

The above relationship avoids the possibility of transverse flow over a grooved no-slip surface, perpendicular 
to an a p plied shear st re ss, wh i ch ha s been analyzed and ob served in a number of studies (see 



Stroock et al 



(j2n02ah : lAidaril tOQ% : IWanel tom-. 



Stroock 



Stroock et al 



et al 



(|2002al ) and 



( 2002b)). Such phenom ena have motivated a tensorial 



Stone et al. 



version of Eq [Tl as discussed in 
length b with a rank-2 tensor b characterizing the surface anisotropy: 



(|2004l ). which replaces the scalar slip- 



up = b • (n • Vu) 



(2) 



where fi is the unit normal (directed into the fluid). iBazant and Vinogradoval (|2008l ) proposed to express 
the tensorial slip condition in the convenient form of a mobility law, where the mean surface normal traction 
T = T • n (for T the Cauchy stress tensor) and some mobility tensor M are used instead of the velocity gradient 
and tensorial slip-length, i.e. 

u'* = M • T (3) 

and discussed general physical constraints on M for different types of fluids and surfaces. 

The work herein analyzes and quantifles this proposed relationship for the case of Stokes flow over a broad 
class of weakly textured surfaces. We focus on horizontally sheared fluid over surfaces with arbitrary periodic 
height fluctuations, and continue the analysis later (Section[8|) to surfaces that also have non-uniform hydropho- 
bicity. To rigorously evaluate the properties of u*, our approach is to construct the flow from a family of 
analytical solutions to the equations of motion, which are superposed as needed to satisfy the order-by-order 
boundary conditions on a topographically complex surface. Section [2] defines the problem and proves the valid- 
ity of Eq [3] by deriving the relationship directly from the Stokes equations for a no-slip surface with arbitrary 
periodic height fluctuations. We then proceed in Section |3] with the crucial and natural question of how to relate 
the details of the surface topography to the mobility tensor. This is considered for the case of arbitrary, small, 
periodic surface corrugations, where a second-order approximate formula for M is constructed from a family 
of Fourier series solutions to the Stokes equations and carrying out a domain-level perturbation analysis up to 
second-order. This result is then used to derive/compute a number of consequential results: analytical results 
in the case of grooved surfaces (Section [S|, flow optimization over surfaces of fixed heterogeneity (Section [HI, 
and statistical results for random surfaces (Section [7|). The limitations of our approximation are also deduced 
and quantified through an in-depth error analysis in Section U] that carries over to the appendix. 



2. Problem setup and basics 

Consider a rigid, periodic surface with height z = H{x,y), with period 2Lx in the x direction, and period 
2Ly in the y direction. Let z — correspond to the bottom of the surface, so that H{x, y) > 0. Above the 
surface is a layer of fluid (of viscosity 77) satisfying the no-slip boundary condition along the surface. It is sheared 
from above, at z — >■ cx), by a horizontal shear traction t — {t^, Ty, 0), which induces a steady flow u(a;, y, z) with 
pressure p{x,y,z). As a convention in this paper, we interchangeably represent horizontal vectors with 2 or 3 
components depending on context — for example r can also represent {Tx,Ty) in planar operations. 

Before continuing any further, let us non-dimcnsionalize the problem to scale rj out of the analysis. Let T 
and £ be arbitrary, flxed units of time and length respectively. Let the unit of stress be rj/T. With these units. 
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all system variables and fields are hereby redefined to be their dimensionless counterparts, e.g. 



u ^ u T/C , T ^tT/tj , p^pT/f] , X ^ x/r , H H/C (4) 

We assume that the Reynolds number is sufficiently small that the fiow satisfies the three-dimensional Stokes 
equations. Under our non-dimensionalization, this gives 

V^u = Vp (5a) 
V • u = (5b) 



and the dimensionless boundary conditions are 



u{x,y,H{x,y)) = 



du 

dz 



(6) 



Far above the patterned surface, the flow must asymptote to a simple linear flow with uniform constant 
pressure. A major goal of this paper is to determine the effective slip — that is, the horizontal vector in the 
asymptotic form 

U(z — > oo) = U'" + TZ (7) 

where 

S = z-{H{x,y)) (8) 

measures the distance above the space-average height of the fluctuations. By definition, a perfectly flat, no-slip 
surface has = regardless of r. When the surface has height fluctuations, is usually non-zero and has a 
direction and magnitude depending on the stress vector r and the surface shape. 

As a consequence of the linearity of the Stokes equations, the effective slip and the traction r must be 
related linearly through a 2 x 2 matrix relationship of the form: 

u'* = M • T (9) 

We refer to the matrix M — M(_ff) as the mobility tensor for the surface H{x,y). The mobility tensor can 
be seen as characterizing the dependence of the net flow properties on the direction and magnitude of the 
applied stress. The linearity of the relationship between and t also means that surfaces with more than two 
symmetry directions must have isotropic mobility. 

It is straightforward to prove EqEl For some fixed surface H{x,y) consider two particular flows: The first 
flow is induced by applying a unit shear traction from above pointed in the x direction, giving rise to some flow 
profile Ui(x, y, z) with corresponding effective slip Ui''. The second flow is induced by a unit traction in the y 
direction, generating a flow profile \J2{x,y, z) with slip U2^. Exploiting the linearity of the Stokes equations, 
the following superposition is an exact solution for the flow u(a;, y, z) induced by an arbitrary shear traction 
T — {ti, T2) applied from above: 

u(x, y, z) = Ti\5x{x,y,z) + r2U2(a;, y, z) (10) 
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Figure 1: (a) Setup for the problem. A fluid is horizontally sheared with stress t at a height of S = oo above a surface with 
arbitrary periodic height fluctuations H{x,y). (b) The induced flow proflle becomes horizontal and linear in z at large heights, but 
is more complicated closer to the surface. Our interest is to determine the effective slip u'', pictured above, which corresponds to 
the extra velocity one would obtain when extrapolating the linear portion down to the mean surface height 2 = 0. 

Consequently, the effective slip u" arising from the shear traction r is 

u-^ = riUi^- + raUa'^ = [Ui'^jUa'] • x (11) 
The matrix [Ui*|U2'*] is om' mobility matrix M. 

3. Computing the mobility tensor 

An exact formula for M in terms of H{x, y) seems very complicated to obtain, and unlikely to have a 
tractable form. In the following sections, we derive an approximate formula for the M tensor in the case of 
small height fluctuations, i.e. 

H{x,y) = eh{x,y) (12) 

for e some small, dimensionless number. 

It is difficult to ascribe a single physical description to e at the outset, since there are an infinite number of 
possible size scales that could be extracted from an arbitrary periodic surface. Which dimensionless size must 
be small for the approximation to be valid? To approach this question systematically, our strategy is to obtain 
the approximate solution for M using domain perturbation theory treating e as a small parameter, and then 
perform an in-depth error analysis to determine how the error of the approximation depends on h{x,y). Once 
the error is quantified in this fashion, it can be used to single out a particular definition of e that is "ideal" in 
the sense that the relative approximation error depends solely on e and no other surface properties. 

We now summarize the major result of this paper. Let h{m, n) compose the set of Fourier coefficients 
corresponding to h{x,y) 

h{x, y)^J2 n)e'^'''-'=+''"y') (13) 

m.n 

where km = rm:/ and fc„ — mr/ Ly are corresponding wavenumbers. Then the effective slip velocity obeys 
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(14) 



where the function M(/i) is defined by 

/ 

M{h) = 



E 



n[m, n) 



A- h 



m ' 



iL2 



/i(m, n) 



E 

(m,n)5^0 

E 

(m,n)7^0 



k k 



/l(r71, 7l) 



2fc2 . 2 

h(m,n) 



(15) 



Hence, — e^M is a second-order approximation to the true mobility matrix M. 

First, observe that the mobihty is a n O(e^) effect for heigh t fluctuations of size 0(e). A second-order leading 
perturbation term was also observed in 



Stroock et al. 



(j2002ar ) for the effective slip of pressure-driven Poiseuille 
flow along sinusoidally grooved walls. Furthermore, the leading-order mobility matrix is symmetric regardless of 
h{x,y) (within the broad range of validity of th e analysis, quantified below). This result was suggested, though 
not proved, in 



Bazant and Vinogradova 



( 20081 ). where the notion of symmetric mobility was supported by a 
statistical diffusion argument and treated as an example of the commonly used Onsager-Casimir relations for 
near-equilibrium linear response. The symmetry of M requires that its eigenvectors are orthogonal, and likewise 
a general surface H should have two orthogonal directions along which the shear stress and slip direction align. 



Stroock et al. 


(2002a 


) and 


Stone et al. 


(2004) 



transverse to the direction of the stress traction, as in 

It is also important to note that with the definitions used here, Eqs [M] and [15] give a negative-definite 
mobility. That is, the effective slip should point away from the traction vector (see Figure [TJb)), which reflects 
the notion that surface fluctuations cause flow resistance compared to a flat surface of the same mean height. 
This outcome can also be seen as dire c tly re lated to our choice of placing the origin for z at the mean height of 
the pattern. In studies such as IWangl (|2003l ). the origin is placed at the peak of the height fluctuations, which 
redefines the slip vector in a way that ensures x • u** > 0, and equivalently redefines the mobility tensor to be 
positive-definite. Positive mobility is desirable from an intuitive standpoint b ecause it more closely represents an 
object slipping on a dissipative, "passive surface" (jBazant and Vinogradoval . 12008 ) . However, for our purposes, 
we find it computationally easier t o perform our analyses with the mean surface height selected as the origin. 



Stroock et al. 



(|2002at ) have also found benefit in setting the origin at the 



Other effective slip studies such as 
surface mean. The remainder of this section details the derivation of Eqs [TJ] and 1151 

3.1. Perturbation expansion 



Our approach utilizes the method of domain perturbations (see iHinchI (|l99ll )). iMiksis and David (jl994l ) 
have also used this technique to study surface flows over small-fluctuation surfaces. Their focus was on two- 
dimensional flows over a surface H{x), and the bulk flow behavior was approximated using asymptotic matching. 
Here, we study the three dimensional problem with bulk behavior solved analytically, and carry out the domain 
perturbation analysis to higher order as necessary to reveal the tensorial properties of the effective slip at the 
surface. 

To begin, we let the velocity and pressure be represented by the perturbation series 



u(x,y, z) = Uo(a;,y, z) eui(a;, y, z) + e^U2(x,y, z) + 0{e') 
P{x,y,z) ^po{x,y,z) + epi{x,y,z) + e^p2{x,y,z) + O(e^). 



(16a) 
(16b) 
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Inserting these into the Stokes equations and grouping same-order terms, it is immediately clear that the 
Stokes equations must be upheld at each individual order. Using the method of domain perturbations, the 
expression of the boundary condition along the periodic surface becomes much easier to deal with by expanding 
in a Taylor series about z — 0: 



u(x, y, eh{x, y)) = = Uo{x, y, 0) + eh{x, y) 



dVLQ 



dz 



+ e ui(x,2/, 0) + eh{x,y) 



z=0 



{eh{x,y)f d'^xiQ 
2 dz^ 

-t-e2u2(x,y,0) + O(e3). 



z=0 



This expression is equivalent to three separate boundary conditions for the first three orders: 



Uo(a;, y, 0) = 
ui(a^,y,0) = -h{x,y) 

"2(2^, y,0) = -h{x,y) 



duo 

dz 

dz 



z=0 



2 = 



h(x,y)'^ a^uo 



dz"^ 



(17a) 
(17b) 

(17c) 



2=0 



The traction condition at — >■ (X) is first-order in magnitude, which implies the following three boundary 
conditions: 



duo 

dz 



dui 

dz 







9U2 

dz 



0. 



(18) 



3.2. Solving for horizontally periodic Stokes flow 

By inspection, we find that the zeroth-order solution is standard simple shear: 

Uo{x, y,z)^Tz po{x,y,z) = Kq. 



(19) 



where Kq is a constant. To compute each of the higher-order terms, we need an analytic general solution to 
the Stokes equations that satisfies the z ^ 00 boundary condition. From such a general solution, we can fit the 
more complicated z = boundary conditions. Due to the horizontal periodicity of the surface, it follows that 
any flow solution be similarly periodic in the x and y directions: 



m,n 

p(x, y, z) = ^ &(m, n, z)e^('="^+*=''^l 



(20) 
(21) 



Substituting these forms into the Stokes equations, a somewhat lengthy general solution for a and b can be found 
(see Appendix [Appendix A I , which ultimately depends on three sets of undetermined coefficients. Whenever 
TO and n are non-zero, a and b both decay exponentially in z as expected, which guarantees satisfaction of the 
upper boundary condition. 
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3.3. First-order term 

To solve for the first-order flow, the undetermined coefficients are solved term-by-term by enforcing the 
bottom boundary condition (Eq ll7b() . which now reads 



ui{x,y,0) = -h{x,y)T. 



(22) 



The result is 



(m,n)#0 \ V^rn + '^n J 



Vi 



Pi{x,y,z) 



ix,y,z)^~hiO,0)ry+ Mm, n)e- ~ rj e^('="-+'="^ 



(m,n)7^0 



(m,n):^0 



(23) 

(24) 
(25) 
(26) 



where Ui — (ui, ui, wi) and Ki is an arbitrary constant. The constant terms in the ui and vi expansions grow 
proportionally to the average surface height ft.(0,0), and provide a first-order isotropic slip contribution, as all 
other terms decay exponentially as z — > oo. 



3.4. Second-order term 

To determine the second-order solution, we take the general solution from Appendix |Appendix A| and enforce 
the z = boundary condition, which now reads 



U2(a;,?/,0) = -~h{x,y) 



3ui 



z=0 



where 



dui 
dz 

dvi 
Ih 

dwi 
dz 



z=0 



z=0 



/p-qrp- n{m,n)e 

(m,n)#0 V '^m "1" '^n 



= ^ z(fc„r,-f A:„r,)/i(m,n)e'('='""+'="'') 

^=0 (m,n)740 



(27) 

(28) 
(29) 
(30) 



as follows from Eqs[ll[21 and US] 

Since we are stopping the perturbation analysis beyond O(e^), there is no need to compute the full solution 
for U2 — one need only determine what the constant term is in its Fourier series solution, Eg 1201 This term will 
be the only one that does not vanish as z — >■ 00, and represents the effective slip. Note that the constant term 
in the Fourier expansion of U2(a;,j/, z) is necessarily equal to the horizontal average (U2) at any fixed z, since 
all other terms average to over the period. Hence, we can compute this constant slip term by averaging over 
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the lower boundary condition, 



U2' = (U2(x, y, 0)) = - ( h{x, y) 



9ui 

dz 



(31) 



Since the term on the right-hand side is the average of the product of two periodic functions of equal periodicity, 
we obtain U2'' from the product of the two Fourier expansions. For any coefficients a{m,n), recall the general 
rule 



h(x^ y) a{m, n)h{m, n)i 



a{m, n)h{m, n)h{—m, — r 



a{m, n) h{m,n) 



(32) 



Applying this rule to the given boundary conditions, we uncover the tensor M.{h) from Eg 1151 giving 

U2'' = -M{h) ■ T (33) 

Combining the constant terms from the first- and second-order flows with the zeroth-order shear flow, we obtain 
the limiting behavior 

u(z ^ 00) = TZ - e h{0, 0) T - e^M{h) ■ t + O(e^) (34) 

= {hix,y)) 

By switching to the variable z, the first-order contribution cancels yielding 



u(z ^ 00) = Tz - e^M{h) ■ T + 0(6^) 
and the slip formula fEq immediately follows. 



(35) 



4. Error analysis 

Ji^.l. The need to compute error hounds 

The significance of the slip approximation (Eqll4|) is tied directly to the error of truncating the perturbation 
expansion at second order. It is of critical importance that we quantify this error to correctly interpret the 
mobility matrix approximation M(e/i) ss — e^M(ft,). 

An example is instructive here. Consider a surface of shallow, square grooves — that is, H{x,y) = €h{x,y) 
where 

, ^ I 1 if 2fc < X < 2/c + 1 for integer k , ^ 

nix, y] — < (36) 
[0 if 2fc + 1< a; < 2fc + 2 for integer fc. 

To use the mobility matrix approximation for this surface (Eq I15p , one computes the Fourier coefficients of 
the above and finds that \h{m,n)\ decays like |l/m|. Likewise, the summand in each of the diagonal entries of 
M must decay as and consequently Mn and M22 are both infinite. This prediction for the mobility is 

clearly erroneous, as it suggests th e effect i ve slip velocity should be infinite (regardless of the choice of e > 0), 
which contradicts known results in 



Wana (|2003l ) and basic intuition. 
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The reason for this failure is because this particular surface introduces an infinite amount of approximation 
error. In view of Ea lT^ the error in u'* must grow as 0{e^), which means the error can be expressed as e^C{h, t) 
for some unknown function C. However, as the example has just demonstrated, certain surfaces h cause C{h, t) 
to become large or infinite, which overwhelms the fact that the prefactor is small. Hence, in the upcoming 
subsection, we bound |C(/i, t)| by some known form |x| Err(/i), so that 

,3 



('M{h)-Tj < e^|r| Err(/i) (37) 

Once the function Err is determined, we can affirm that the mobility matrix approximation M « — e^M is valid 
for any surface where the error bound is small compared to the approximate slip, i.e. 

^ , e^Err(/i)|T| e^ETT(h) , , 

Relative Error Bound = ^' ' « < 1 (38) 

I - e2]V[(/i) . x| t'^\yi{h)\ 

Such a relation gives an a priori method for ruling out surface shapes for which the slip matrix approximation 
should fail. 

Below, we show that psp can also be used to determine an ideal definition of e in terms of properties of the 
surface. Up to this point, we have described e merely as a "small dimensionless number", and h{x,y) is defined 
as H{x,y)/e. A specific formula for e as a function of the surface H(x,y) is physically desirable, wherein the 
size of e immediately relates to the accuracy of the approximation. 

^.2. Bounding the approximation error 

The second-order fiow solution is a superposition of analytical solutions to the Stokes equations. Hence, the 
PDEs are always satisfied exactly. Likewise, the error of the slip approximation arises solely from error in fitting 
the boundary condition u(x, y, e/i(a;, y)) = 0. 

Appendix [Appendix B| constructs a bound on this error, which leads to the following bound on the error of 
the second-order slip approximation: 

- (-e'M(/i) • r) I < e^n /im|x|(/im |VV;^|,, + \Vh\l,) (39) 

where we use the subscript m to denote a function's maximum value over all x and y. The constant k, is 
dimensionless and independent of r and the surface shape; a first calculation gives k < 55. Based on the 
analysis in Appendices [Appendix B| and [Appendix C[ we speculate it is not possible to bound the error in a 
general fashion using fewer than the first two derivatives of h. While this decrees a qualitative level of tightness, 
it is clear from the appendices that quantitatively tighter bounds could be written using non-local norms and 
surface integrals. We prefer the above, as it is expressible using basic quanti ties. Recent work on bounding the 



error of Reynold's lubrication approximation for thin-channel Stokes fiows (jWilkeningj . |2009|) has found some 
similar dependences on the norms of the derivatives of the surface. 

Equation IMl reveals that the error is unbounded for any surface where |VV/i|^,j or |V/i|Af oo. This 
result gives the following crucial limitation: The approximation M(e/i) w —e^'M.{h) should not be used 
for surfaces with corners or vertical slopes. This point, which might not be surprising to those familiar 
with domain perturbation techniques, provides an explanation for why the square grooves example of the last 
subsection failed. The mobility matrix formula can only be reliably applied to functions h{x,y). Note that 
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the error does not depend at all on the higher-order derivatives of h (order > 3), and as such the approximation 
can be accurate even if the higher derivatives of the surface shape are discontinuous. 

Equation 1551 can also be used to read off a formula for the function Err that was sought after in Eg 1571 



Err(/i) = K /iM (ft-M |VV/i| 



M 



(40) 



Observing EqsHOlandlTSl we note that the relative error bound (Eg can be expressed solely in terms of the 
actual surface H{x, y), and completely independent of e, since all powers of e can be brought into the arguments 
of the functions M and Err, and eliminated by using H = eh: 



Relative Error Bound 



e3Err(/i) _ Err(iJ) 
e2|M(/i)| " \M{H)\ 



(41) 



This last expression is a dimensionless ratio determined entirely by the given surface H{x,y). While it may 
not be simple to interpret physically for a general surface, its direct connection to the error suggests an ideal 
formula for choosing e given H{x,y) would be 



Err(i/) 

mm 



kHm{Hm I VVffU, + |Vff II,) 



E 

(m,n)#0 

E 



fc2 



H(m, n) 


CO 




2 


H{m, n) 





(42) 



H(ra, n) 

\(m,n)7^0 

The relative approximation error can never exceed eideai, hence all surfaces with a small value for eideai 
are well-described by the mobility formula Eq [151 If we let e — eideai in Eq[37l and consequently h{x,y) = 
-ff (cc, y)/ eideai, then the slip relation can be rewritten as 



-e^deaMh) ■ T 



< e: 



deal 



mh)\ \t\ 



(43) 



The definition of eideai reduces to a simple, expected form in canonical cases. For example, consider the 
surface 

H{x,y) = asin(a;/&) (44) 

for constants a and b. We expect that the relevant dimensionless quantity in this case should be a/6, in the 
sense that the slip approximation over this surface should improve as the height to wavelength ratio decreases 
to zero. Evaluating (j42|, the relevant small quantity according to the error bound is eideai = (8^/^/5) (a/6), 
which is proportional to a/6 as expected. 

To conclude this discussion on error analysis, we emphasize the major points. Most importantly, the approx- 
imate slip given by Eqs 1141 and 1151 can be inaccurate for certain surfaces. Approximation error grows with the 
maximal curvature and maximal slope of the surface. To determine the accuracy of the approximation, there 
are two mathematically identical representations that have been discussed (Eqs [35] and . The representation 
shown in Eq [Ml is most applicable when the user wishes to decide the value of e based on a loose observation 
of the characteristic surface height. Such an approach is useful when computing the mobility properties of a 
broad family of surfaces, all with similar characteristic height e. To survey the range of surfaces, one need only 
vary the scaled surface h, taking care to ensure the range includes only those h shapes that keep Eq [Ml within 
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one's error tolerance. On the other hand we also present the error bound in a second form, that of Eg 1431 which 
is in terms of a specific definition of e extracted directly from the surface H{x, y). This error representation is 
best when there is only one surface H to be studied, since Cideai is in the truest sense the small parameter of 
interest, which immediately determines the relative accuracy of the approximation. However, if a comparison 
over many surfaces is desired, the need to recompute Cideai for each surface H can complicate the analysis. 



5. Mobility over surfaces varying in only one direction 

Our main result above (Eq llSp is useful to reveal general properties of the effective slip and, within the 
error bounds, to allow its computation for an arbitrary surface. However, the various series for the elements of 
M can be cumbersome to evaluate. Of course, it is easy to evaluate them for surfaces with simple sinusoidal 
perturbations, which have Fourier series truncated after a small number of terms, but this offers little analytical 
insight into the dependence of the mobility tensor on the shape of the surface. To more clearly see this 
dependence, this section focuses on a simpler family of surfaces. 

Observing Eq 1151 we note that the formula reduces significantly for the class of surfaces varying in only one 
direction. Letting the variation direction be x, these surfaces appear as parallel groove patterns with shape h{x) 
and mobility 



M w -e^^A(h) = -e^P ( " " 1 , for /3 = 2 V fc^ ' ~ 



=1 



/i(m) (45) 



This shows that the slip length, which measures distance from the mean surface height to the equivalent no-slip 
plane, is generally twice larger for perpendicular versus parallel shearing. For anisotro pic Stokes flow, factors o f 



two like this are not surprising — it is reminiscent of similar results for striped pipes (jLauga and Stond . l2003l ) 



and the clas sical result that a rod sediments twice as fast in creeping flow if aligned vertically, rather than 



horizontally (|Batcheloij . ll97Cl ). 

Let us evaluate M in closed form for a nontrivial subset of grooved surfaces, whose Fourier series have an 
infinite number of terms, and vary continuously in shape from sinusoidal to sharply peaked. Consider surface 
shapes h{x) of the type 

h{x) = ^'^^ ~ , cj,{x) = (46) 

(Pmax - (Pmin 1 + + 2a COS X 

and the unit of length is £ = Lx/n for simplicity, so that km = m. Regardless of the parameter < a < 1, 
h{x) has a maximum of 1, and a minimum of 0. As shown in Figure^ a controls the shape of the surface. For 
a <C 1, the surface height is a small perturbation from a single- mode sinusoid. As a — > 1, the surface develops 
wide deep valleys around the minimum height at a: = 2mT separated by tall, narrow peaks at x = (2n + l)7r. 

The Fourier coefficients of h{x) can be obtained directly from those of (f>{x). Set z ~ to obtain a rational 
function f{z), separate into partial fractions, and expand in geometric series (since \az\ — \a/z\ < 1) to obtain 
the Laurent series of f{z), which equals the desired Fourier series on the unit circle: 
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2n 4k (tn 



Figure 2: The surface shape h{x), defined in Eg 1461 under different values of the parameter a. 



Consequently, 



By Eq[ 



h{m ^ 0) = \ ( "^'"J' = ^(-a)l™l (47) 

m=l ^ ^ ^ ^ T?x=0 

We observe that for this specific family of surfaces, /3 does not vary with a. However, the observed slip 
properties will vary, since the average surface height, at which the slip calculation applies, does change with 
a. To understand the slip behavior more clearly, we use Eq to compute a second-order approximation 
to the location of the equivalent no-slip plane (denoted z^ "^', i.e. the height of the plane having zero mean 
velocity) corresponding to each a under both parallel and transverse shear. The mean surface height obeys 
e{h) = e(l — a)/2, and consequently, for any e and unit applied shear stress, we have to second-order 

- + ^ (49a) 

1-a 



,N.S. 



This qualitatively states that as a increases, the effective no-slip plane descends (regardless of the shear direc- 
tion). In physical terms, the fluid content in the groove pattern increases with a, which increases the lubrication 
of the surface. The directionality of the surface comes in at second-order, where the added flow resistance of 
shearing transverse to the grooves is apparent. 

As a increases, it is important to keep in mind the accuracy of the approximation. The second derivative 
of h and the square of the first derivative both diverge as (1 — a)^^ as a — >• 1. Hence, Eg [39l implies the error 
bound on z^'^' must also diverge as a — > 1, 

Azf^-, Az^/-<e3^^^0(l) (50) 

Such a bound is useful for detecting qualitative behavior of the accuracy, however, as is common with asymptotic 
methods like domain perturbation, the actual error is commonly well under the bound. 

To test this, we compute a near-exact flow solution for parallel shear over the family of surfaces (Eg H^. 
Since the velocity is always parallel to the grooves, the problem reduces to solving the Laplace equation for the 
in-line component v{x^z) under the boundary conditions Vz{z = oo) = t = 1 and v{x^eh{x)) = 0. This can be 
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10" 



1-a 



Figure 3: The true error of the second-order approximation for z^-^- relative to the size of the second-order correction (i.e. 
vs. the value of a. The rate of divergence near a = 1 is notably lower than our conservative error bound oc (1 — a)~^. 



solved as a superposition of separable solutions 

oo 

v{x, z) = z + Am cos{mx) exp(— mz) (51) 

■m=0 

where the Am are chosen to satisfy the no-slip boundary condition. To compute a numerical solution, we solve 
for the first 1000 terms Am by enforcing the boundary condition on 1000 equally spaced x values from to tt. 

Flows are solved using a values ranging from 0.15 to 0.99 in increments of 0.03. For each a, the flow is 
computed using seven e values decreasing from 10^^ to 10^^ logarithmically. Our goal is to determine how 
accurate Eg I49al is compared to the finest term of the approximation, the second-order correction Hence, 
for each a and e pair, we compute the relative second-order error — the difference between the near-exact 
numerical solution for z^"^- and the second-order approximation (Eg I49a|) . scaled by the size of the second- 
order correction term As expected from the error bound fEa [50]) . the relative second-order error comes 
out as being proportional to e and diverges to cx) as a — > 1. However, the proportionality constant is in general 
much smaller, and diverges more slowly than the bound. The results are plotted in Figure [31 From the plot, 
we see that if e = 0.1, for example, then as long as a is less than about 0.8, the second-order prediction for the 
no-slip plane is accurate to within one tenth of the size of the second-order correction term. 



6. Optimizing slip properties 



With a direct approximation for the mobility matrix and bounds for the error, we now attempt to answer 
the question: What surface shapes minimize/maximize the effective slip? Problems of this type, relating to flow 
optimization with rigid interfaces, have been looked at p rimarily in three dimensio ns, through elegant analyses 
of optimally porous structures for fluid permeability (see I Jung and Torguatd (|2005r )l. Here we study essentially 
a reduced dimensional problem of optimizing the effective planar slip with respect to a rigid boundary. We study 
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(a) 



(b) 



Figure 4: (a) The periodic surface of fixed variance that restricts flow the least (i.e. maximal forward mobility) has long wavelength 
sinusoidal grooves oriented along the direction of the stress, (b) The surface of fixed variance that restricts flow the most (i.e. 
minimal forward mobility) has sinusoidal grooves of the shortest possible wavelength oriented across the direction of the stress. 



the case of macroscopic shearing, thou gh in the opposite hmit of a thin-channel flow, similar wall-optimization 
studies have recently been performed (jFeuillebois et all 120091 ) . 
Let us define ^ ^ 

Forward mobility = " ' j^j^^'^^'' ~ - ' M(/i) • (52) 



which measures how much slip occurs in the direction of r per unit shear stress. As is evident from Eg 1151 the 
most forward mobility occurs when the surface is perfectly flat: If h{x,y) = const, then M = and = (as 
expected from a flat no-slip surface). For all other surfaces, M is necessarily positive definite and the forward 
mobility must be negative. In simple terms, fluctuating surfaces always resist flow more than flat surfaces. 

We next investigate the related question of how to maximize/minimize the mobility given that h{x,y) must 
have some fixed level of heterogeneity. We quantify the heterogeneity through the variance of h{x, y): 



Yaiih) = {ihix,y)-{hix,y))f)= ^ \him,n) 



6.1. Maximal forward mobility 

Suppose e is fixed, = Ly = L, and h{x, y) is constrained to have a fixed variance . We now derive the 
surface shape h{x,y) that maximizes the mobility. This is essentially a Lagrange multiplier problem where the 
unknowns are the coefficients \h{m,n)\. Without loss of generality, presume t is aligned with the x direction. 
Then our goal is to minimize 



X • M • X = Mil = ^ 



2kl 



(m,n)7^0 



h{m,n) 
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subject to the constraint 



The critical points occur only for 



Var(/i) = |ft,(m,n) 



(53) 



V|.,_,|MiiW = AV|,.^,_,l Var(M 
for Lagrange multiplier A. Evaluating this form gives 



2fc„ + fc„ 



2 I jU2 



h{m,n) — 2X h{m,n) for all (m,n) 7^0 



(54) 



Hence, any surface constituting a critical point of the Lagrange multiplier problem can only contain modes with 
the same value of (2fc^ + fc^)/-\/fc^7+~^- ^'^y critical point containing a certain mode 6*^'^°^+'^''''' must, by 
Eqs[53]and[15l have 



Mil 



+ H 2 



(55) 



Of these critical Mn values, the minimum value is Afn = ira'^/L, which occurs only for {a, (3) — (0, ±1). 
Together with the fact that h{m,n) = h{—m, —n)*, we obtain 



e't h{0, -1) 



'''I', and h{m^O,n^±l) = 



for some constant phase (/). The minimizing surface given our constraints is then 



h{x,y) 



72 



V2 



aV2 cos 



(56) 



Thus, at fixed variance, the least resistance to forward flow occurs when the surface is a single-mode sinusoid 
and the fluid is pushed in the direction along the grooves. Moreover, the ideal wavelength of the surface is the 
longest one allowable for the given periodicity. 

6.2. Minimal forward mobility 

Before determining the surface shape with the least forward mobility, an important caveat must be included: 
To wit, we constrict the allowable bandwidth of the surface. For some positive K, consider only scaled surfaces 
of the form 

h{x,y)= ^ /i(m, n)e^('^-^+'="2') (57) 

Capping the bandwidth is the simplest way to avoid surfaces that violate the second-order accuracy of the 
perturbation expansion. Naive observation of Eg 1551 would suggest that minimal mobility occurs as a,/? — > 00, 
as this maximizes the critical values of Mn. But, in view of the error bound Eg 1391 we see that such a solution 
violates the accuracy requirements, as it corresponds to a surface with infinite frequency and hence unbounded 
slope and curvature. By considering only surfaces of a finite bandwidth, we guarantee there exists a single, fixed 
e so that all critical points of the Lagrange multiplier problem correspond to sufficiently accurate solutions. 
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This feature is not a major setback in the study of reahstic surfaces, where material structure usually possesses 
some inherent "grain size" placing natural limits on the sharpness of corners. 

This being said, we reduce to the problem of minimizing the forward mobility over the set of surface shapes 
h with Var(/i) = ct^, and maximal wavemimber K. Let 

T=IKL/ti\. (58) 

Then by inspection of Eq [55l the largest critical Mn value is 2Ttt<t'^/L, which occurs only for /? — 0, and 
a = ±r. The corresponding surface is 

Hence, the surface with the least forward mobility has sinusoidal grooves oriented perpendicular to the direction 
of the shear stress. The wavelength of the grooves should be the shortest possible, given the fixed surface 
periodicity and fixed bandwidth. In light of the result from the previous subsection, a very natural symmetry 
arises in the solutions to maximum/minimum mobility. The overall directionality of the grooves in each case 
is not surprising — the transverse and parallel dis tributions are ubiquitous in the optimization of physical 



properties for a variety of materials (jTorauatol . l2002l ) — though the groove shape is perhaps less obvious. 



7. Random height fluctuations 

Based on the prior section, an appropriate follow up is to determine the mobility properties of a (2L-periodic) 
random surface under an identical set of constraints, those being: 

1. The surface is real, and therefore /i(m,n) = h{—m,—n)*. 

2. The surface has Var(/i) = a^. 

3. The bandwidth of h is finite, with maximal wavenumber K. 

Let A be the set of integer pairs (m, n) with the property that + fc^ < K , m > when n > 0, and 

m > whenever 71 < 0. Under these definitions, any surface obeying constraints (1) and (3) above must be of 
the form 

hix, y) ^ h{0, 0) + J2 ('*("^' n)e*('='"^+'="^) + h{m, n)* g-'C^'-^^+'^-y)) (60) 

{m,n)eA 

and the variance constraint (2) is 

J2 2\h{m,7^)\^^a\ (61) 

(m,n)GA 

Observing the slip matrix formula Eq 1141 the mobilities depend only on the magnitude of the Fourier 
coefficients. This suggests a simple, broad probability distribution for sampling the surfaces as shall be described 
next. First, index the members of A by pi through p|yi| and define jj = 2|/i(pj)p. The region of space in the 
variables {71, ...,7|^|} that fulfills the constraint on the variance is expressed by 
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m 



Figure 5: The integer grid points within an upward semicircle, excluding the non-positive m-axis, compose the set A. Up to an 
additive constant, any periodic surface of finite bandwidth can be uniquely described by determining h{m, n) on A. 



< 71 < cr^ 
< 72 < cr^ 



71 



< 73 < -71-72 



(62a) 
(62b) 
(62c) 



\A\-2 

< 7|A|-i <<y^ - ^3 
j=i 

7|A| = cr^ - ■ 



(62d) 
(62e) 



Let us denote by 57 the higher dimensional tetrahedron described by the equahties above on 71, • • • ,7m_i. 
Each surface h obeying our constraints is equivalently one point in the region fi. So, one probabihty distribution 
we can write is 



Prob of /i € with 7's between 7^ and 
7j + ^7j fo'' "^^ch 1 < j < 1^1 — 1 




n 



(63) 



j=i 



Under this distribution, the probabihty of choosing some surface within a subset of occurs with probabihty 
proportional to the volume of that subset. Using an overbar to denote the weighted average by this distribution, 
basic integration gives the identity 
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J = l 3 

for any function /. Applying this to EqllJl one finds 



\A\ 2 1-41 

E/(p.)^.- - ^E/(p.) (64) 



/-^ /),2 I 1,2 /i-2 4_ h.2 



(65) 



Next, each of these sums over A must be computed. Each sum in the matrix above can be approximated by 
the semicircular integral 



pKL/7T />7r 

f{m,n)= / / f {r cos 9, r sin 9) rdrdO 



f{m,n) ^ ^ 

(m,n)eA 

due to symmetry along the m-axis and the fact the /(0,0) = for each summand. Then we may write 

Lit 

2lT 



M=rTT ^-Sfs^cr^^f-r-rrl- (66) 

|A| I ^ 2tt\A\ ^ ^ 



where 1 is the 2x2 identity tensor. Lastly, the number of elements in A can be computed with a semicircular 
area approximation (excluding the missing point at the origin): 

\A\ - l7:iKL/nr 1 

giving the final result 

where the last approximation holds for KL large enough. 

As expected, a random surface of fixed variance and bandwidth has isotropic slip properties. Moreover, the 
forward mobility of a random surface is approximately —e^a^K which fits nicely into the ordering we find for 
the most and least optimal surfaces of fixed variance as discussed in section [SJ 

Forward mobility: -e^a^^ < ^_eV^ < -2€'^a'^^^ 

^ ^ ^ Random 
Most 

Or, perhaps the ordering is easier to see in terms of the relative sizes of these three quantities, which is closely 
approximated by 

Most : Random : Least = - (1 : F : 2r) 

But we reiterate that the result for mean mobility over a random surface depends entirely on the choice of the 
probability distribution. We have selected a particular one, which is simple enough to perform the necessary 
calculations, and which appears relatively unbiased toward any particular subset of surfaces. 
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8. Surfaces with fluctuating scalar Navier slip 



Up to this point, tlie focus has been on surface patterns due solely to height fluctuations. Now, let us 
expand the analysis to surfaces that have fluctuating scalar slip properties, such as composite surfaces where 
each material has different hydrophobicity properties. 

Such surfaces can be described with a scalar Navier slip boundary condition, which relates the fluid slip 
along the surface to the fluid shear-rate at the surface 

n^^b{x,y)^. (68) 

A no-slip surface corresponds to a Navier slip coefficient of 6 = 0, whereas a traction-free surface corresponds 
to = oo. In this section, we analyze the macroscopic flow properties of surfaces with periodic slip coefficient 
eh{x,y) measuring small variations from no-slip, as well as possible height fluctuations eh{x,y). Applying the 
Navier slip condition along the surface eh{x, y) gives the boundary condition 

n{x,y,eh{x,y))^eh{x,y) .{"t^' ^j, ■ Vu ^ ^ (69) 



{-eK, -ehy, 1 



z — eh(x^y) 



Adopting the perturbation series representations of u and p from section 13.11 one expands in e to yield the 
following term-by-term conditions for the flow at z = 0: 



uo(x,y,0) = 

^i{x,y,0) = {b{x,y) ~ h{x,y)) 



9uo 

dz 



"2(2;, y,0) = -h{x,y) 



dvLi 

dz 



z=0 

h{x,yY d'^VLQ 



dz^ 



+ 6(x, y) [ h{x, y) + ( "^^ ^ 



z=0 

5uo 



z=0 



\ — 

^ dy 



z=0 



1 • Vui 



z=0 



The z — 00 conditions are the same as those in section [3TT1 The term Uo is solved, as before, by 

^o{x,y,z) = Tz. 

Substituting this result into the above boundary conditions for the other orders, one finds 



Ui{x,y,Q) = {b{x,y) ~ h{x,y)) 
U2{x, y, 0) = {b{x, y) - h{x, y)) 



dvLi 

dz 



(70) 
(71) 

(72) 

(73) 

(74) 

(75) 
(76) 



Comparing to sections l3.3l and l3.4l it is now evident that these boundary conditions are the same as the previous, 
except with the switch h{x, y) h{x, y) — b{x, y). Hence, the solution up to second-order can be inferred directly 
by modifying the results of the previous sections. 

Keeping with the convention that z = at the bottom of the surface, we apply the switch and obtain the 
result 

u(z -> 00) = TZ - e{h{x, y) - b{x, y))T - e^m{h - b) ■ t + 0{e^) (77) 
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Writing, as before, in terms of the variable z = z — e{h{x,y)), one obtains the sUp relation 

= e{b{x, y))T - e^M.{h - b) ■ t + 0{e^). (78) 

Up to second-order, the effective slip of a surface with small scalar slip and height fluctuations is the same as 
that of a no-slip surface except with /i — > /i — 6 in the argument of M, plus an isotropic first-order term due to 
the average of b. 

8.1. A simple optimization for surfaces with height and Navier slip fluctuations 

Consider a surface with height and slip fluctuations of the same periodicity. The forward mobility, as 
computed from Eq[78l is 

^e{b{x,y)-e' \ ' . 79 

Holding the average of b{x,y) fixed, the mobility can be maximized by observing the second-order term above. 
As per the definition of M (Eq llSp . the mobility matrix is positive definite if b(x,y) — h{x,y) ^ const, and 
otherwise. Consequently, the mobility is maximized when b{x,y) ~ h{x,y) — const. For any given h{x,y), the 
optimal choice of y) is 

y) - {b{x, y)} = h{x, y) - {h{x, y)} (80) 

which gives a mobility of e{b{x,y)). Physically speaking, this means the forward slip is maximized when the 
peaks of the surface h arc more hydrophobic, and the valleys of the surface are more hydrophilic. 



9. Conclusion 

This work has derived a second-order accurate formula (Eqs ll4l and ll5p describing an effective local boundary 
condition for shear flows over small-fluctuation periodic surfaces. The formula represents a tensorial mobility 
law, and is easily extended to include surfaces of both non-uniform hydrophobicity and height changes f Eg [751) . 
We have gone to great lengths to quantify the error of the approximation, so as to provide firm guidelines for its 
appropriate usage. Within these guidelines, the formula was optimized in a Lagrange multiplier framework to 
derive the no-slip surfaces of fixed variance that max;imize/minimize the forward mobility. These extremal cases 
were favorably compared with the mobility of a random, fixed-variance surface as computed by summing over 
some unbiased distribution. We have also performed a simple optimization that instructs the optimal coupling 
between height and scalar slip fluctuations when both surface effects can take place. 

In the future, we hope to augment our analysis to include the possibility of a surface charge profile. This 
addition would be useful in electro-kinetic applications to help understand and predict electro-osmotic fluid 
transport. We also continue from a mathematical perspective to work on enhanced solution methods that may 
improve the breadth of applicability of the mobility formula. 



Appendix A. Appendix: General solution to the Stokes Equations 

In this appendix, a general solution to the Stokes equations is given for use in solving the first- and second- 
order terms in the perturbation expansion displayed in Section 13.11 It is expressed as a Fourier series in the 
horizontal dimensions with three undetermined coefficient sets: 
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u{x,y, z) = a + e 



(7n,n)^0 



jU2 I 1.2 



B{m,n) C{m,n) 



U2 I 1.2 



J(kr„x+k„y) ' 



3\/fcm + kn 



'ikn (1 + 2zyjkl^ + fc2^ A(m, n) 



kmkn{\/ k'^-^ + fc2 2z(fc2 + fc2 3A;2 •\/fc2 + fc^ 2/s^z + 2/0^ (-^/fc^ + fc2 fc^ 2:) 

B(m, n)H — — C(m,n) 



h'^ J_ (^2 



fe2 I 



(m, 71)7^0 



3 + 2z^/k^~\-~kfA A(m, n) — 2ikmzB(x, y) — 2zfc„zC(TO, n) 



p{x,y,z) = S+ ^ e'(fe™^+fc.y) 

(m, 71)7^0 



'i'v/^m + k'^A{m, n) — 4ik„iB{m, n) — AiknC{m, n) 



This system can be inverted to uniquely determine the constants a, /?, 7 and the coefhcient sets A{m, n), B{m, n), C{m, n) 
in terms of the chosen lower boundary condition on the velocity, u{x,y,0). The general solution automatically 
upholds the known boundary conditions for the first- and second-order perturbation terms: 



u{x,y,z) = u(x + L^,y + Ly,z) , p{x,y, z) = p{x + L^,y + Ly, z) , 



du 



0. 



Appendix B. Appendix: Determining error bounds 

This appendix derives the error bound displayed as Eq in Section 14.21 To begin, define the following 
norm on scalar, vector, or tensor fields f (x, y) over a 2D periodic domain: 



lf(^,2/)ll. 



1 



4:LxLy J-L^ 



|f(x,?/)|2 dxdy 



(B.l) 



We reserve | • | for the absolute value (if applied to a scalar) or Euclidean norm (if applied to a vector or matrix) . 
The above is directly proportional to the compact support L2 norm on functions. The choice of this norm 
greatly simplifies the analysis since it it connects directly to an inner product, and consequently an array of 
useful theorems. 

The exact flow solution in our problem satisfies the no-slip condition along the surface, i.e. u(a;, y, eh{x, y)) = 
0. To measure how much our second-order flow approximation differs from the true solution along the surface, 
we use the norm to measure the surface error: 



Surface Error = \\uo{x,y, eh{x,y)) + eui{x, y,eh{x,y)) + e'^U2{x,y,eh{x,y))\\^^y 



(B.2) 
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From Taylor's Remainder Theorem, we know 



Ui (a;, y, eh) = Ui {x, y, 0) + eh 



U2{x,y,eh) = U2{x,y,0) + 



dui 



\h(x,y) 



=0 JO 



{eh{x,y) - z)-Q^ dz 



dU2 

dz 



dz. 



Inserting these forms into the above, and using the order-by-order boundary conditions to cancel terms, one 
obtains 



Surface Error 



< e 



{eh{x,y) - dz 



{eh{x, y) — z) „ ^ dz 



dz^ 



+e' 



dz 



dz 







dz 



dz 



(B.3) 



by the triangle ininequality. 

To place an upper bound on the error, we must calculate tight bounds on each of these two integrals given 
any surface h(x, y). This is a highly nontrivial task. As is derived in Appendix [Appendix C[ with judicious use 
of Parseval's Theorem and the Cauchy-Schwartz inequality, we can show that 



{eh{x,y) — z) „ „ dz 



dz^ 



<e'hljK,\T\ ||VVMa;,y)L,, 



iih{x,y) 



9U2 

dz 



dz 



<ehMK2\T\{\Vh\M \\'^h{x,y)\\^,y + hM || VV/i(a:, 



(B.4) 
(B.5) 



x,y 



where Ki and K2 are order one, known, dimensionless constants independent of the choice of e, r, or h{x,y), 
and for any function /(cc, y) we adopt the notation 



Jm =maxf{x,y). (B.6) 
The bounds on i?^^' and i?*^^' applied to Eq IB. 31 give 

Surface Error < k e^hM\T\{hM \\\"\/h{x,y)\\x^y + |V/i|Af \\'Vh{x,y)\\x,y) 

<Ke^hM\T\{hM |VV/i|j,,f + |V/(,||f) EE Surface Error Bound (B.7) 

where we define k = Ki -\- K2, which is another known order-one constant. 

Next, we make the connection between the surface error and the error in the predicted slip. The exact 
flow solution must have the velocity go to zero along the surface. The second-order approximation that we 
have found has some fallacious extra velocity along the surface that arises due to truncation error, and the 
surface error is a space-average measurement of this erroneous surface speed. If we let u(x, y, z) represent the 
true solution, u°'PP{x,y,z) = Uq(x, y,z) + eui{x,y, z) + e^U2(a;,y, z) be our second-order approximation, and 



22 



A(a;, ?/, z) = u{x, y, z) — u°p^'(x, y, z), it follows that 



\\A.{x,y,eh{x,y))\\^ y < Surface Error Bound (B.8) 

Since both u(a;, y, z) and u°'Pp{x, y, z) obey the Stokes equations exactly, the difference A.{x, y, z) must also 
obey the Stokes equations. Moreover, since both u and u^^'^ satisfy the same traction boundary condition as 
z — >■ oo, then the difference flow A(x,?/,z) must asymptote to a constant uniform flow (zero shear traction) for 
large z. The speed in the uniform flow region is precisely the error in the approximation for the effective slip. 
Our goal is hence reduced to determining the maximum possible flow speed that can occur at large z in a Stokes 
flow A(a;, y, z) that obeys Eg IB.8I and zero shear tractions at z = oo. 

While the total surface error is described by ()B.7|) . the distribution of this extra velocity over the surface 
is undetermined, providing the only degree of freedom in this maximization. We state without proof that the 
maximal large-z uniform flow arises for 

A.{x,y,z) — (Surface Error Bound) e (B.9) 



where e is any horizontal unit vector. In other words, A has a large-z uniform flow of maximal speed when the 
extra surface velocity represented by the surface error is equally distributed and uniformly directed along the 
surface, so that A{x, y, z) is globally uniform. This choice of A saturates Ea IB.8I 

Recalling the definition in (|B.7|) . the above selection of A{x,y, z) implies the error in the second-order slip 
approximation is subject to 



which is our desired result. 



Appendix C. Appendix: Bounding the integrals R^^^ and R^'^^ 



The purpose of this appendix is to rigorously prove the bounds on the error quantities and shown 
as Eqs IB. 41 and IB.5I in Appendix [Appendix B[ 

Appendix C.l. The error 

We begin with the bound on R^^\ First, we observe by the Cauchy- Schwartz inequality that 



< y^€h{x,y)]^ J {eh{x,y)~z) 

< \/ehM\ / (e/iA/ - z) 



feh(x,y) 



{eh{x,y) - z) 



dz^ 



ih{x,y) 



eh{x,y) 



12 dz 



dz^ 



dz 



dz^ 



dz 
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We insert this result into the definition of R^^\ and expand the function norm: 



Ve/i^W J (e/iM - z) 



dz 



\ 



{ehM - zf- 



92 Ui 



/e/i 



M 



Jo 



92 Ui 



dz^ 



dz. 



(C.l) 



The final line follows from rearranging the order of integration. In order to apply the next step, we must first 
compute \d'^ui/dz'^\'^, which is straightforward because ui{x, y,z) is completely known. The result is lengthy, 
but can be written somewhat compactly as 



a2ui 



a^2 



a^2 



+ 



9^2 



+ 



9^2 



J = {1,2,3} \(m,n)#0 / 

where for each choice of m and n, q_m,n and Pm,n are the vectors defined by 

Ctm,n — (1) Z\/k^ + A;2 , zkjn, zkn) and Pm,n = {kmt ^n' V^femfcrn kjn\/k^ + fc2, kn'\/k^ + fc2) (^.2) 

and for each j, A(j-| and B(j') are constant, dimcnsionlcss, 4x5 matrices of order one in size. 

Now for the crucial step: Apply Parseval's Theorem to convert the function norm to a discrete sum, 



52ui 


2 i-L^ pLy 


a2ui 


dz"^ 




dz"^ 



dxdy 



E 

J={1,2,3} 



L 



^ ^ ' (Qm,n * ^(j) ' Pm,n ; Qm,7l * ^{j) ' Pm,n) 

» \(m,n)7tO 



dxdy 



i={l,2,3} (m,n)#0 

< E E M'|q-,nl'(|A(i)|' + |B(,-)|')|p„,„|2e-2-v^^^|Mm,n)|' 

j={l,2,3} (m,n)7tO 

The last line follows from the Cauchy-Schwartz inequality. To express the result more simply, define 



"= E |Ao-)l' + |Bo-)| 

j={l,2,3} 
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Then, expanding the norms on the vectors p and q, wc have 

2 



9^2 



x,y 



< a|rp a + ^H^kl + ^kl)){2kf^ + 2kt + Aklkl)e-'^V^^. h{m,n) 



(C.3) 



With a bound on the integrand in (|C.1[) now estabUshed, we move on to evaluate the integrah 



ehM 



{ehM - zf 



dz 



x,y 



< 



{ehM-zfa\T\' J2 il + z'i2kl + 2kl))i2kt,, + 2kf, + iklkl)e-'^V^^\h{m,n) 



dz 





2 




h{m, n) 




/ 






Jo 



(ehM - z)2e-2V^-^+'=?. dz 



ehM 



+2{kl + kl) / [ehM - zfz^e-^^V^^^ dz 



(C.4) 



Now we bound the two z integrals appearing in the last expression. One can show by basic calculus that for 
any 2: > 0, 

e-2-\/^OT < 1 and z2g-2.yfeIT^ < 



e^{kl, + kiy 



(C.5) 



Hence, 



and 



ehM 



{ehM - z)2e-2V=™+'=' dz< [ehM - zf dz 



ehM 



3u3 



ehM 



{ehM - z)'z'e-^^VWl dz < 



'Hkl + kl)J, 



ehi. 



{ehM — zY dz 



4^3/, 3 



3e2(A^ + fc2)- 



Substituting these bounds in Eq (|C.4|) . we obtain 



ehM 



{ehM - z) 



9z2 



dz < e^/iL|x|2A'2 ^ {kt + kt + 2klkl) h{m, n) 



(C.6) 



25 



where we define Ki = ■\/2a(l + 8/e'^)/3. Now, observe that 



ikt+kt, + 2klkl)\h{m,n) 



(m,ji)7tO (m,n)5^0 {m,.n)^0 



(m,n)7iO 



I — fc^/i(m, ri) +2 I — fcmA:„/i(m, n) 

(m,n)^0 (m,n)7^0 
2 



(m,n)5<tO 



(m,n)7tO 



(m,n)5^0 

2 



d^h{x,y) 


2 

+ 


d'^h{x,y) 


2 

+ 2 


d'^h{x,y) 






9y2 




dxdy 



(C.7) 



Parseval's Theorem is invoked in the penultimate hne to swap each discrete sum for the function norm of a 
Fourier series. Substituting Eq lC.7l into Eq lC.61 we may rewrite Eg IC. II as 



M 
2,2 



d^h{x, y) 


2 

+ 


d'^h{x,y) 


2 

+ 2 


d'^h{x, y) 


2 


dx"^ 




dy"^ 




dxdy 





e'hijK^\T\\\V^h{x,y)\\ 



giving us our final result. 



Appendix C.2. The error 

To compute the second error term i?^^)^ one must have the full solution for the second-order velocity 
U2(x,y, z). Let C„, C^,, and C.^, represent the non-constant-term Fourier coefficients for U2(x,?/, 0), obtain- 
able via convolution. That is, 



U2{x,y,0) = -h{x,y) 
V2{x,y,0) = -h{x,y) 
W2{x,y,0) = -h{x,y) 



dui 
dz 

dvi 
dz 

dwi 
dz 



z=0 



z=0 



-X • M • T -I- C„(m, n)e 



(m,n)7^:0 



(C.8) 
(C.9) 
(C.IO) 
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Now, inserting the boundary conditions into the general solution form, one obtains: 



(m,n)#0 V +K J 



V2 



W2 



^ e-"V^OT ((^A:2^ + fc2C^ - imCu - inC,)z + C^) e'('='""+'=-^') (C.13) 



(m,n)7^0 



The derivation of the bound on i?'^^) follows a similar sequence as the proof in the last subsection for R^^'^ 
Steps will be abbreviated accordingly. First, the Cauchy-Schwartz inequality allows us to bound the integral: 



dz 



dz 



M 1 



9U2 



dz 



dz 



As in the previous section, the next step is to insert this result into the definition of i?^^) ^nd expand the 
function norm, ultimately giving us 



< A eh 



M 



ehn 



92 U2 

dz 



dz. 



(C.14) 



Taking the z derivative of U2 above, we find that 



dU2 



dz 



^ ^ I ^ ^ ^ni,n ' {^ni,n ' -^(j) * Pm,n 7 Qm,n ' -^(j) ' Pm,n ; 
i={l,2,3} \(m,n)#0 

Qm.n ' (j) ' Pm,n) 



where the vectors Pm,™ and (\ra,n are defined in Eqs lC.21 Cm.n = {Cm Cu, Cui), and the matrices D(j), Eqj, and 
Fq) are constant, order one, 4x5 matrices for each j. 

As was done before, we next apply Parseval's Theorem to this result, giving 



^■2/ j = {l,2,3} (m,n)5^0 



^m,n ■ (Qm,n ' -^(j) ' Pm,n 5 Qm,n ' -^(j) ' Pm,n 1 Qm,n ' {j) ' Pm,n) 
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We apply the Cauchy-Schwartz inequality to this result, and simplify /expand algebraically: 



dU2 



dz 



— ^ ^ |Cm,n| |Pm,n| |Qm,n| 



I 1,2 



Ij. |2|p |2 ^^m + + ^^m^n ^-2z^/k'^„+kl 



(m,n)7iO 



E 2/3 |c™,„np„,„|2(fc2^ _^ ^2)^-2.VeT^ 
(m,ra)5^0 

5] 2/3 |c„,„p(l + 2z2(fc^ + fc^))(fc^ + fc^)^-2.yfeM 



(m,n)7^0 



where 



i={l:2.3} 

The integral in Eg IC.14I can now be bounded: 



'E/,n I + |F|- 





dU2 


Jo 


dz 



dz< J2 2/3|c„,,„p(fc2^ + fc2 

(m,n)/0 



ehM . /"e^ivj , 

e-^^V^^?^dz + 2{kl + kl) I z^e-^-V^^^dz 







where K2 = ■\/2/3(l + S/e^) is obtained by bounding the two integrands per (jC.5|) . 

Substituting this bound into Eq IC.141 apply Parseval's Theorem again, keeping in mind the definition of 
(C„,a,C^) from EqsEH 



(m,n)5^0 



(m, 71)51^0 



= e^M^2 



\ 



ehMK2 



dui 



z=0 



V /i(a;,y) 



9z 



2 = 



x,y 
= ehMK2 



d_ 
dy 



h{x,y) 
{Wh{x,y)) 





dui 
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dz 






dui 

dz 


+ 

z=0 


h{x,y) 



dui 

dz 



z=0 



< ehMK2 \Vh\ 



M 



i9ui 












dz 


z=0 


+ h-M 


dz 


z=0 





(C.15) 



Computing \d\ix/dz\'^ at z = 0, the solution takes the form: 



dvLi 



dz 



z=0 



j = {l,2,3} \(m,n)740 



For 



^771,71 I f) 



and t,„ ,j — (^7711 kn) 
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and Gq) and H(j) constant, dimcnsionless, order one, 3x2 matrices. Parseval's Theorem and the Cauchy- 
Schwartz inequahty then give: 



dui 




dz 


z=0 



2|A |2 

It 



— ^ ^ T I'^l |^m,n| |^m,n 
2 /;2 , ,2 



/i(m, n) 



(m,n)5^0 



27 |r| 



fc„) /i(TO,n^ 
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9/i 




dx 
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dy 





27 kniV/i||^,, 



(C.16) 



for 

By similar means, we can show 

2 



7= E |Go-)l^ + |Ho)| 

J = {1,2,3} 







dz 


z=0 



< E 7 (fc™ + A')kns™,„nt™,„|' /i(m,n) 

(m,n)^0 

= 27 |t|2 ^ I ui , ol2 ,„2 



= 27 |t| 



/i(m, n) 



( 




2 
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2 

+ 2 




') 








9y2 




dxdy 





27 \T\'\\VVh\ 



2 



(C.17) 



Applying (|C.16p and (|C.17|) to ()C.15|) . and absorbing -y/27 into the definition of K2, we obtain the final result: 

< ehMK2\T\ (\Vh\M ||V/i||^_^ + hM II VV/i||^ ,^ 
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